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ABSTRACT 


The oscillation of a gas bubble produced as a result of underwater explosion could 
cause the severe whipping damage on nearby marine vehicle. The effects of rigid boundary 
curvatures to explosion gas bubble oscillation behavior in underwater were investigated. 
The analyses were conducted using a multimaterial Lagrangian-Eulerian finite element 
code, MSC/DYTRAN. The incident shock wave pressure, bubble pulse pressure, gas 
bubble radius and period were calculated for the case of detonation of a charge near the 
curved rigid boundary. The results were compared for the case of free field bubble 


oscillations. 
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I. INTRODUCTION 


When a charge detonates in underwater, it emits energy to the surrounding water and 
converts the original material into a gas sphere (called the bubble) at very high temperature 
and pressure. The emitted energy, which has about half of the energy of the explosion, 
appears as a shock wave. During the transmission of the shock wave to the surrounding 
fluid, the initial high pressure in the gas sphere is considerably decreased, but it is still 
much higher than the surrounding pressure, which is atmospheric plus hydrostatic 
pressure. As the shock wave passes each water particle, the water particle attains a radial 
velocity. This motion of the water particle behind the shock wave front is essentially 
incompressible flow after the pressure in the water drops to near hydrostatic pressure. Thus 
it is reasonable to use incompressible flow theory to describe the bubbles motion and 
behavior. Initially, the bubble has a large outward velocity and its diameter increases 
rapidly. The expansion continues for a relatively long time, and the internal bubble pressure 
decreases gradually, but the motion persists because of the inertia of the outward flowing 
water. At a later time, the pressure in the bubble falls below the equilibrium value 
determined by the surrounding pressure and the pressure defect brings the outward flow to 
stop and the bubble begins to contract at an increasing rate. As bubble approaches 
minimum radius, the bubble produces a pressure pulse, called the bubble pulse. The bubble 


pulse can give a radial velocity to water particle and cause damage just as the initial shock 








wave can. The bubble contains about half of the total explosion energy[Ref. 1]. Explosions 


in air emits a shock wave which has the most of explosion energy, and do not form a 





bubble like an underwater explosion. In underwater explosion analysis, the bubble can be 


as important as the initial shock wave. 


Explosive gas bubble behavior in the free field case and several simple boundary cases 
was studied by Chisum[Ref.2]. The effects of rigid boundary curvatures on gas bubble 
oscillation behavior in underwater explosions was investigated. The auaivnes were 
conducted using the multimaterial Lagrangian-Eulerian finite element code MSC/DYTRAN. 
The incident shock wave pressure, bubble pulse pressure, gas bubble radius and period 
were calculated for the detonation of a charge near various curved rigid boundaries. The 


results were compared with the free field bubble oscillation. 











II. NUMERICAL COMPUTER CODE 


In this study, the materials analyzed are fluids ( gas and seawater ). The Eulerian 
finite element technique, which is frequently used for analyses of fluids or materials that 
undergo very large deformations, and in which the grid points are fixed in space and the 
material moves through the Eulerian mesh (see Figure 1) is used. The mass, momentum 
and energy of the material are transported from element to element. The Eulerian method 
calculates the motion of material through elements of constant volume with respect to 
time(see Figure 2). MSC/DYTRAN is capable of finite element analyses using to 9 Eulerian 
materials. In the MSC/DYTRAN code, the Eulerian mesh is defined in exactly the same 
manner as a Lagrangian mesh. General connectivity is used, so the Eulerian mesh can be of 
an arbitrary shape and have an arbitrary numbering system. This offers considerably more 
flexibility than the logical rectangular meshes used in other Eulerian codes. However, the 
use of an Eulerian mesh is different from that of the Lagrangian type. The most important 
aspect of modeling with the Eulerian technique is that the mesh must be large enough to 
contain the material after ilimeuiies A basic Eulerian mesh acts like a container and, 
unless specifically defined, material cannot leave the mesh. Stress wave reflections and 


pressure built-up can develop from an Eulerian mesh that is too small for the analysis. 








Figure 1. 











Transportation of Mass, Momentum and Energy through Eulerian Elements 
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Figure 2. Motion of Material through Eulerian Elements 








Il. MODELING AND SIMULATION 


In order to investigate the effects of rigid boundary curvature on explosion gas bubble 
behavior, the free field case which has no boundary was analyzed first, to provide basic 
information about the bubble. Cases in which the explosion takes place near a flat ngid 
boundary were then simulated, to compare to cases with curved rigid boundaries. Finally, 
curved rigid boundary cases were analyzed. In this study, the charge modeled is a 10.24kg 
cylinder of TNT, 20cm in a diameter and 20cm height, submerged in a seawater at a depth 
of 1000m. The TNT was modeled using a Jones-Wilkins-Lee (JWL) equation of state, 
with standard equations and parameters for TNT [Ref. 4]. With this equation, the pressure 
in the “burned fraction” of the explosive material is related to the specific internal energy 
and density by : 
oh Ry 

n 


+B(1- rae " +enpE (1) 


aon 
p= Al-—)e 
k, 
where N= /p, 
O )= initial density 
E = specific internal energy (per unit mass) 


A,B, w, R,, R, are constants for the explosive 
The explosive was modeled as a homogeneous mass of TNT. The equation of state 


parameters used were : 

















9, = 1630 kg/ m° 
E = 20076 J/kg 

A =3.712xX10" Pa 
B =3.231x10° Pa 
@ = 0.30 

R, =4.15 

R, =0.95 

d = 6390 m/s 


In order to model the seawater in which the explosion occurs, a polynomial equation of 
state was used. This state of equation relates the pressure in the fluid to the acoustic 


condensation w and the specific internal energy by; 


p=apt+a,p’ +a,u°+(b, +b,+b,u’)p,E (tL >0) 
(2) 
p=a+(b, +b )poE (LL < O) 
where LL = (P-Po)/Po 


0,= initial density 
E = specific internal energy (per unit mass) 


a,,a,,a,,b,,b,, b, are constants for the fluid 
and the upper equation applies to a fluid in a compressed state, while the lower applies to a 


fluid in a expanded state. Constants were determined by Chisum [Ref. 2]. 
A. FREE FIELD CASE 


The free field case was investigated first. As deep submerged charges are known to 


undergo little vertical migration due to gravity [Ref. 1], a quasi one dimensional finite 








element model, in which gravity is neglected entirely, was used for this analysis. The 


problem thus has spherical symmetry. The geometry of this problem is shown tn Figure 3. 


p (atm) = 0.1013 MPa 






p (atm) = 0.1013 MPa 
p= 1025 kg/m? 
y=1000m distance r = 0.71m g = 9.80665 m/s” (1) 


p.. = p(atm) + pgy 
(= 10.15MPa) 


charge (10.24kg TNT) pressure measure point (r : standoff distance) 
Figure 3. | Deep Spherical Bubble Problem Geometry 


The overall geometry of the finite element model for this quarter symmetry 
axisymmetric free field analysis is shown in Figure 4, from a three dimensional 
perspective. Figure 5 shows the complete finite element model for this analysis from a two 
dimensional perspective. The distance of 400m from the charge to the boundary was 
chosen — at this distance, reflected wave form the boundary did not affect the results 
during this analysis. A close up view of the lower night corner of Figure 5, which includes 
the “central” which was a refined region near the charge and “transition” region which 
was used to connect the central and outside regions, and a small portion of the “outside” 
region which can be thought of as “storing” energy and momentum during bubble 


expansion, which then effects bubble collapse after the momentum field reverses direction, 








MN 








Figure 4. Overall Model Geometry for Quarter Symmetry Axisymmetric 


Finite Element Model for Free Field Case 
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Figure 5. Quarter Symmetry Axisymmetric Finite Element Model for Free Field Case 
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Figure 6. 
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Close Up View of Finite Element Model in Area near Charge 


for Free Field Case 
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is shown in Figure 6. The black area represent the charge. For simplicity in model 
generation, the two angled faces in the models which appear to meet at the y-axis are 
actually offset a very small distance (+ 5um) ;they thus make an angle of +1 degree not 
from the xy-plane at z=0, but from the xy-plane at z=+ 5m and z=- Sum, respectively. 
This small offset can not be seen on any reasonable scale (itis small even in comparison 
with the other dimensions of the smallest elements), and permits the model to be meshed 


using only hexahedron elements. 


Figure 7 shows the bubble radius time history. The maximum radius (r,) 1s about 
0.71m, and the period T, is about 14.7 milliseconds. Figure 8 shows the total pressure 
time history. This pressure was measured at 0.71m (same as maximum bubble radius) 
from the charge center in the horizontal direction. The shock wave pressure and first bubble 
pulse peak pressure were about 270 MPa and 46 MPa, and occurred at 0.5 and 14.7 
milliseconds, respectively. The second and third bubble pulse magnitudes are not clear, but 
are still higher than the surrounding pressure (10.15 MPa). Empirical equations for TNT 


[Ref. 1] predict a bubble radius, period and dynamic pressure given by : 





maximum bubble radius 
| ea 
R. =3.382 -) 3 m 
— X Dal o (m) 
bubble period 
j 
w3 
T  =2.064——~ (second) (3) 
(D +10)° 
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Figure 7. Bubble Radius Time History for Free Field Case 
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Figure 8. Pressure Time History for Free Field Case 
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dynamic pressure 


1 


W 3 
Prax = 52.118(——)"" (MPa) 


where D = charge depth (1m) 
L= standoff distance (m) 
W = charge weight (kg) 


These empirical equations give a maximum bubble radius of 0.71 m, a first oscillation 
bubble period of 14.06 milliseconds, and a total shock wave pressure of 205 MPa. These 
values are approximate and not much different from the values predicted by the finite 
element analysis. Figure 9 shows radial velocities of a fluid particle at the same point 
where the pressure was calculated. This computation was done using a spherical wave 
equation ind the pressure time history given in Figure 8, but using dynamic pressure vice a 


total pressure. The spherical wave equation used is shown below [Ref. 5]. 
1 1 ¢t 
u(t) =u,(t)+— [p(t)-po.t+—] paTxT (4) 
Poo Pore 


where iu, (t) = radial velocity of fluid particle 
u,(t) = initial radial velocity of fluid particle 
9, = initial fluid density (1025 kg/m°) 
c, =acoustic velocity in seawater (1500 m/s) 
r =standoff distance 


p = dynamic pressure 


From this equation it is evident that the radial velocity of a fluid particle depends on 


both the pressure at that time and the pressure time history. 
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Fluid Particle Velocity {m/s} 





E+0i 


Time (msec) 


Figure 9. Radial Velocity Time History of Fluid Particle for Free Field Case 
B. FLAT RIGID BOUNDARY CASE 
Figure 10 shows the bubble problem geometry for the flat rigid boundary case. 


p (atm) = 0.1013 MPa 







9 =1025 kg/m° 
g = 9.80665 m/s” (1) 
p.. = p(atm) + pgy 

(= 10.15MPa) 


y=1000m distance r = r, (0.71m) 


a 


h* charge (10.24 kg TNT) pressure measuring point 


( rigid flat boundary ) h* : normalized standoff distance 


Figure 10. Bubble Problem Geometry for Flat Rigid Boundary Case 





The typical overall geometry of a typical finite element model for a flat rigid boundary case, 
and the complete finite element model, are almost same as shown in the figures for the free 
field case; the distance from the charge center to the outer boundary is also the same as in 
the free field case. A view of the finite element model near the charge is shown in Figure 
11; the black area represents the charge. The modeling concept is the same as used in the 
free field case, and was also used for the subsequent cases. In this analysis, the pressure 
measuring point is same as in the free field case. The distance between the charge center 
and the rigid boundary surface is defined as a normalized standoff distance h* = h / Irs 
where h is the standoff distance from the charge center to the rigid boundary surface, and r, 
is the maximum bubble radius for the free field case(0.71 m). Analyses were conducted for 
h* values of 1.16 and 2.0 . Obviously, the free field case is the case for which h* = 


infinity. 


Figure 12 shows the bubble radius time histories for these cases. The maximum 
bubble radii are 0.67 m for both cases. The first bubble periods (T) for the cases where h* 
= 1.16 and h* = 2 are 16.6 and 15.5 milliseconds, respectively. The pressure time histories 
for both cases are shown in Fi gure 13. The shock wave pressures appear at the same time 
as in free field case and the magnitudes are also the same. At 1 and 1.7 milliseconds, 
pressures of 86 and 38 MPa appear in Figure 13 for the cases where h* = 1.16 and 2.0, 
respectively. Using the empirical equation for pressure and Snell’s law to calculate the 
pressure at the pressure measuring point via bottom reflection, with a constant acoustic 
velocity in seawater of 1500 m/s, the reflected total pressures are computed as 75.4 and 
46.6 MPa for the both cases, and occur at 1.2 and 1.96 milliseconds, respectively. The 


secondary pressure peaks at 1 and 1.7 milliseconds are therefore bottom reflected 
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Figure 11. Close Up View of Finite Element Model in Area near Charge 


for Flat Rigid Boundary Case 
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(a) case for h* = 1.16 
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(b) case for h* = 2.00 


Figure 12. Bubble Radius Time History for Flat Rigid Boundary Case 
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(b) case for h* = 2.00 


Figure 13. Pressure Time History for Flat Rigid Boundary Case 
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pressures. Bubble pulses are 31 and 38 MPa at 16.6 and 15.5 milliseconds. As h* 
increases, the maximum bubble radius increases and the bubble period and reflected 
pressure decrease, approaching the results from the free field case. If standoff distance 
h>c,T,/2, where c, is the acoustics wave velocity (1500 m/s) and T, (14.7 
milliseconds) is the bubble period for the free field case, the flat ngid boundary cases will 
be treated as the free field case. So in cases where " > 15.5 are considered as the free field 
case. This is because the bubble behavior will not be affected in the first bubble period by 


the boundary when h 2 15.5. 


Figure 14 and 15 show the maximum bubble radius, non-dimentionalized dividing by 
r, (radius for free field case ), and the first bubble period, non-dimentionalized dividing by 
T, ( period for free field case ), versus the normalized standoff distance h*. The cases for 


h*=1 and 4 were calculated by Chisum [Ref.2]. In these figures, h* = 30 indicates h* = 


infinity. 
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Bubble Radius Ratio ( 
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Figure 14. Bubble Radius Ratio vs. Normalized Standoff Distance h* 
for Flat Rigid Boundary Case 
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Figure 15. Period Ratio vs. Normalized Standoff Distance h* 
for Flat Rigid Boundary Case 
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Figure 16. Bubble Problem Geometry for Rigid Boundary Case 


with Positive Curvature 
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C. CURVED RIGID BOUNDARY CASE 
1. Cases With Positive Curvature 


Figure 16 shows the bubble problem geometry fora ri gid boundary case with positive 
curvature. The finite element model for an axisymmetric rigid boundary with positive 
curvature case is shown Figure 17. The distance from the charge center to the outer 
boundary is 400m. Figure 18 shows the “central” region of the model; black area 
represents the charge. The empty semi-circular arc represents the boundary surface with a 
normalized arc radius of r* = R // r, where R is the actual arc radius. Analyses were 
conducted for in which r* is equal to 1, 2, 4 and 10, and the normalized standoff distance 


h* is 2. 


The bubble radius time histories for each case were shown in Figure 19. Table 1 lists 
the maximum bubble radius and bubble period. Figure 20 shows the pressure time 
histories. Table 2 indicates the shock wave pressure, reflected pressure and first bubble 
pulse. In these analyses, the bubble acts like the free field case when r* is less than 2. As 


r* increases, the bubble behaves like the flat rigid boundary case with the same h*. 


Figure 21 and 22 are rearranged to show the maximum bubble radius and first period, 
non-dimentionalized by dividing by r, ( radius for free field case ) and dividing by T, 
(period for free field case), respectively, versus the normalized arc radius r*. In these 


figures, r* = 30 indicates r* = infinity. 
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Figure 17. Axisymmetric Finite Element Model for Rigid Boundary Case 


with Positive Curvature 
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Close Up View of Finite Element Model in Area near Charge 


Figure 18. 


for Rigid Boundary Case with Positive Curvature 
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(b) case for r* = 2, h*¥ = 2 
Figure 19-1. Bubble Radius Time History for Rigid Boundary Case 
with Positive Curvature 
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(c) case for r* = 4, h* = 2 
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| (d) case for r* = 10, h*¥ = 2 
Figure 19-2. Bubble Radius Time History for Rigid Boundary Case 
with Positive Curvature 
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| (b) case for r* = 2, h* =2 
Figure 20-1. Pressure Time History for Rigid Boundary Case 
with Positive Curvature 
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(d) case for r* = 10, h*¥ =2 
Figure 20-2. Pressure Time History for Rigid Boundary Case 
with Positive Curvature 
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Table 1. Max. Bubble Radius and First Bubble Period for Positive Curvature 


| | max. bubble radius (m) | first bubble period (msec) 


















shock wave reflected bubble pulse note 
pressure (MPa) | pressure (MPa) (MPa) 
| none | free field case 





flat boundary 
(h* = 2) 
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Fi gure 21. Bubble Radius Ratio vs. Normalized Arc Radius r* for Case 


with Positive Curvature 


Period Ratio (=T/Tf) 
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Figure 22. Period Ratio vs. Normalized Arc Radius r* for Case with Positive Curvature 








2. Cases With Negative Curvature 


Figure 23 shows the bubble problem geometry for a rigid boundary case with negative 
curvature. The finite element model for an axisymmetric rigid boudary with negative 
curvature is shown Figure 24. In this case, the charge is enclosed by the curved ngid 
boundary and the distance from the charge center to the boundary is not 400m. Figure 25 
shows the “central” region of the model; black area represents the charge. Analyses were 


conducted for values of r* is equal to 4, 7, 9 and 13, and a normalized standoff distance h* 


is 2. 


The bubble radius time histories for each case are shown in Figure 26. Table 3 lists the 
maximum bubble radius and first bubble period for each case. Figure 27 shows the 
pressure time histories. Table 4 indicates the shock wave pressure, reflected pressure and 


first bubble pulse. 


In these analyses, the charge was enclosed by a rigid boundary and distance from the 
charge to the boundary was very small compared to the other cases in which the distance 
was 400m. ( note that the case when r* = O is impossible ) For the negative curvature 
boundary cases , the results are affected by multiple reflection paths. This is already seen in 
the graph for the r* = 4 case where there is a large amount of “noise”. In the r* =4 case, the 
boundary is so close that it significantly effects bubble growth and collapse. As r* 
increases, the bubble radius and the bubble period are seen to increase, approaching the flat | 


boundary case results for large r* values, as expected. 
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Figure 23. Bubble Problem Geometry for Rigid Boundary Case 


with Negative Curvature 


Table 3. Max. Bubble Radius and First Bubble Period for Negative Curvature 
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Figure 24. Axisymmetric Finite Element Model for Rigid Boundary Case 


with Negative Curvature 
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Close Up View of Finite Element Model in Area near Charge 


Figure 25. 





for Rigid Boundary Case with Negative Curvature 
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(a) case for r* = 4, h* =2 


Bubble Radius (cm) 
E+01 





0 0.5 1 1.5E+01 
Time (msec) 


(b) case for r* = 7,h* =2 
Figure 26-1. Bubble Radius Time History for Rigid Boundary Case 
with Negative Curvature 
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(c) case for r* = 9, h* = 2 
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(d) case for r* = 13, h*¥ =2 
Figure 26-2. Bubble Radius Time History for Rigid Boundary Case 
with Negative Curvature 
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(a) case for r* = 4, h*¥ =2 
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(b) case for r* = 7, h* =2 
Figure 27-1. Pressure Time History for Rigid Boundary Case 
with Negative Curvature 
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(c) case for r* = 9, h* =2 
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(d) case for r* = 13, h*¥ = 2 
Figure 27-2. Pressure Time History for Rigid —s Case 
with Negative Curvature 

















Table 4 Shock Wave, Reflected Pressure and Bubble Pulse for Negative Curvature 


shock wave reflected bubble pulse 
pressure (MPa) | pressure (MPa) (MPa) 
free field case 
(reference) 





flat boundary 
(h* = 2) 





Figures 28 and 29 show the non-dimentionalized bubble radius and period as a 


function of the normalized arc radius r*. In these figures, h* = 30 indicates h* = infinity. 


D. BUBBLE SHAPE 


The shapes of the bubble for the free field case, the flat rigid boundary cases with h* = 
1.16 and 2 and the case with positive boundary curvature with r* = 10 and h* = 2 are 
plotted in Appendix A at times t of 0.0, 0.50, 0.80, 0.90 0.95, 0.98, 0.99, and 1.00T, 
where T is the first bubble period for each case. The free field case and the rigid boundary 
cases with h* = 2 have similar shapes in which the bubble remains nearly spherical up to 


1.00T. The case with h* = 1.16 shows a heart shaped bubble after time 0.90T. The center 


38 























of mass of the bubble in this case is also seen to move closer to the boundary at late times. 


These related phenomena are due to the effect of the boundary on the fluid flow 


surrounding the charge. 
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Figure 28. Bubble Radius Ratio vs. Normalized Arc Radius r* 


for Case with Negative Curvature 
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Figure 29. Period Ratio vs. Normalized Arc Radius r* for Case with Negative Curvature 





IV. CONCLUSION 
A. FREE FIELD CASE 


The free field case gives the basic information about the bubble behavior in an 
underwater explosion. The empirical equations are valid for this case, and the spherical 


wave equation can be used to calculated the fluid particle velocity with reasonable accuracy. 
B. FLAT RIGID BOUNDARY CASE 


The flat rigid boundary case is a special case of a curved rigid boundary case with the 
normalized arc radius r * equal to infinity. The normalized standoff distance h* affects the 
bubble behavior significantly. As h* increases, the bubble behavior approaches that of the 
free field case. If the standoff distance h> c,T,/2, there is no boundary effect. When h* 
is small enough, the bubble is attracted to the boundary surface and the bubble becomes 


heart shaped during late stages of the bubble collapse. 
C. CURVED RIGID BOUNDARY CASE: 
1. Positive Curvature 


In this case, the normalized arc radius r* influences the bubble behavior. When r* 
is small enough, the bubble behavior is similar to the free field case. And as 1* increases, 


the bubble behaves as in the flat rigid boundary case. 


Al 








2. Negative Curvature 


The normalized arc radius r* has a similar effect upon the behavior of the bubble 
as in the positive curvature case. However when 1* is small, there are multiple reflection 


paths which effect the bubble and make the analysis of this case complicated. 
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V. RECOMMENDATION 


This paper describes the effects of rigid boundary curvature on the bubble behavior in 





an underwater explosion. Since the boundary is rigid, this does not represent an actual 


case. Next time, the flexible boundary case needs to be investigated. 
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APPENDIX : BUBBLE SHAPES 
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Figure 30-1. 
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Figure 30-2. 
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(d) at Time = 0.90T 


Bubble Shape for Flat Rigid Boundary Case with h* = 1.16 
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Figure 3 1-2. 
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(d) at Time = 1.00T 


Bubble Shape for Flat Rigid Boundary Case with h* = 1.16 


Mitts 




















$-—-> X one x 
(a) at Time = 0.0T (b) at Time = 0.50T 
2 ' 
$c ¥ ee ‘f 
(c) at Time = 0.80T (d) at Time = 0.90T 


Figure 32-1. Bubble Shape for Flat Rigid Boundary Case with h* = 2 
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(a) at Time = 0.95T (b) at Time = 0.98T 





(c) at Time = 0.99T | (d) at Time = 1.00T 
Figure 32-2. Bubble Shape for Flat Rigid Boundary Case with h* = 2 
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(a) at Time = 0.0T 


(c) at Time = 0.80T 


Figure 33-1. 





(b) at Time =0.50T 





(d) at Time = 0.90T 


Bubble Shape for Rigid Boundary Case 


for Positive Curvature with r* = 10 and h* = 2 
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